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Abstract 

Elastic models of the glass transition relate the relaxation dynamics and the elastic properties of structural glasses. 
They are based on the assumption that the relaxation dynamics occurs through activated events in the energy landscape 
whose energy scale is set by the elasticity of the material. Here we investigate whether such elastic models describe the 
relaxation dynamics of systems of particles interacting via a purely repulsive harmonic potential, focusing on a volume 
fraction and temperature range that is characterized by entropy-driven water-like density anomalies. We do find clear 
correlations between relaxation time and diffusivity on the one hand, and plateau shear modulus and Debye-Waller 
factor on the other, thus supporting the validity of elastic models of the glass transition. However, we also show that 
the plateau shear modulus is not related to the features of the underlying energy landscape of the system, at variance 
with recent results for power-law potentials. This challenges the common potential energy landscape interpretation 
of elastic models. 
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1. Introduction 

The relaxation dynamics of supercooled liquids oc¬ 
curs through local particle rearrangements whose en¬ 
ergy cost is related to the elastic properties of the mate¬ 
rial. This suggests the existence of correlations between 
the elasticity and the dynamics of supercooled liquids. 
Indeed, at the local scale particle mobilities have been 
found to be related to local elastic constants ll]]. At 
the macroscopic level these correlations have stimulated 
the formulation of elastic models of the glass transition, 
which relate the relaxation dynamics and the elastic¬ 
ity of glass formers. For instance, according to Dyre’s 
shoving model ||2l, the relaxation occurs through local 
volume-preserving events that allow the system to tran¬ 
sit from one potential energy minimum to a different 
one. If one estimates the energy barrier separating two 
energy minima within a parabolic approximation, this 
approach leads to the relation 

T oc e.x^(GpVJk^f 

between the relaxation time and the plateau shear mod¬ 
ulus, Gp, where Vat is an atomic volume element that is 
assumed to be temperature independent, T is tempera¬ 
ture and is Boltzmann’s constant. See Refs, yll^lst] 
for a discussion regarding the interpretation of Gp. Re¬ 
cent numerical results on polymer melts ||3l and on sys¬ 
tems of particles interacting via inverse power law po¬ 
tentials 1^, showed the possibility of connecting the 
plateau shear modulus to features of the energy land¬ 
scape of the system Si. Indeed, these studies found 
Gp to be related to the fluctuations of the inherent shear 
stress, Gp - G'® = Here is the 

shear stress measured after quenching a system to its 
inherent structure, i.e. the nearest potential energy min¬ 
imum. Together, Dyre’s shoving model and the results 


of Refs. iSIS lead to a relation between the relaxation 
time on the one hand, and features of the inherent en¬ 
ergy landscape on the other. 

A related approach to connecting the dynamics to 
elastic properties via features of the energy landscape 
can be motivated by an analogy with the Lindemann 
melting criterion for periodic crystal structures. Here 
a relaxation event is considered to occur through local 
rearrangements that take place when the mean squared 
vibrational amplitude of a particle (u^), which is its 
Debye-Waller (DW) factor, crosses some threshold a^, 
where a is of the order of the particle size. If this pro¬ 
cess requires an energy barrier AE oc k^Ta^Ku^) to be 
overcome, one recovers the Hall-Wolynes equation 

T OC exp(a^/2{M^)); 

this connects the structural relaxation time to a short- 
time elastic property, the DW factor. This approach has 
recently been 1^ generalized by introducing a proba¬ 
bility distribution for a^, and successfully tested against 
experimental and numerical data, including both strong 
and fragile glass-formers. The main message from this 
work is that there is a one-to-one correspondence be¬ 
tween the relaxation dynamics and the DW factor. 

These two approaches for understanding glassy re¬ 
laxation times that we have described above are close 
related to each other because the DW factor is hxed 
mainly by the shear modulus of the material E3l- We 
will therefore refer to them collectively as “elastic mod¬ 
els”. 

In this paper we consider the applicability of elas¬ 
tic models to suspensions of particles interacting via a 
harmonic potential. Similar hnite range purely repul¬ 
sive potentials are of interest as model for the interac¬ 
tion of macroscopic particles such as bubbles, foams 
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and microgels, whose dynamics exhibit glassy features 
at high concentration and/or low temperature. In addi¬ 
tion, these potentials are also of interest for being able 
to give rise to water-like density anomalies at high den¬ 


sities 111 lLll2Lll3Lll4ll5n . in spite of their manifest sim¬ 


plicity. The possible applicability of elastic models to 
these systems is of particular interest because elastic 
models are based on an energy landscape interpretation 
of the dynamics, while density anomalies have been ra¬ 
tionalized by entropic arguments ll2ll . 

We will show that in finite range repulsive systems 
the shear modulus estimated from the properties of the 
underlying energy landscape overestimates the plateau 
shear modulus, which implies that the energy landscape 
properties are poorly correlated with the elasticity. Nev¬ 
ertheless, we do find correlations between relaxation 
time, diffusivity, plateau shear modulus and DW fac¬ 
tor that are consistent with those predicted by the elastic 
model of the glass transition. Our results indicates that 
while elastic models correctly describe the relaxation 
dynamics of the system, their interpretation in terms of 
features of the energy landscape should be reconsidered. 


2. Model 


We consider a poly disperse mixture of N = 10^ har¬ 
monic disks of mass m, in two dimensions. Diameters 
are uniformly distributed in the range [Dmin : 
with the difference Dmax - Dmm between the largest 
and smallest diameter being 82% of the mean diameter 
(Dmax + Dinin)/2 SO that the distribution is fairly broad; 
this is necessary to prevent crystallisation at high vol¬ 
ume fractions. Two particles i and j with average diam¬ 
eter D - {Di + Dj)l2 and at a distance r interact via a 
potential 


v(r) = 



if r < D 
if r > D 


( 1 ) 



Figure 1: (color online) Mean square displacement for a system of 
N = 1000 particles at T = 0.18 and 0 = 1, 1.15 and 1.3. 


where Va, , Faij and rdj are the a-components of the ve¬ 
locity of the i-th particle, of the force between particles 
i and j and of the separation between them. The tran¬ 
sient shear modulus is related to the decay of the shear 
stress fluctuations, G{t) - j^(cr^j.(0)cr„(f)). We ex¬ 
plore the features of the underlying energy landscape 
by repeatedly minimizing the energy of the system via 
the conjugate-gradient protocol to find the instanta¬ 
neous inherent structure, and measuring the inherent 
shear stress cr^^. The latter is computed via Eq. 
where all velocities are set to zero. The transient in¬ 
herent structure shear modulus is defined as = 

^(cr^y(0)cr^®(r)), where T is the temperature of the 
parent liquid; G*®(0) is then the same as G^^ defined 
above. 


3. Density anomalies 


In the following, lengths, masses and energies are ex¬ 
pressed in units of Dmax, ni and of e, respectively, and 
the density is expressed via the volume (or rather, area) 
fraction cj) - N{A)IL?. Here L is the system size, (A) 
the average particle area, and N the number of parti¬ 
cles. We have performed molecular dynamics simu¬ 
lations d at fixed volume, temperature and particle 
number, integrating the equations of motion using the 
Verlet algorithm, and constraining the temperature via a 
Nose-Hoover thermostat. The shear stress is defined as 
the off-diagonal term of the stress tensor. 


N 


- - 




'xijFyij 


\ i=l 


i*J 


( 2 ) 


Systems of harmonic and Hertzian particles are char¬ 
acterized by water-like density anomalies nil [Hill 
IHIIl , including a non-monotonic variation of the dif¬ 
fusivity upon isothermal compression, as well as a neg¬ 
ative thermal expansion coefficient. As an example, we 
show in Fig. [T] the mean square displacement for three 
different values of the volume fraction, (p - 1, (p - 1.15 
and (p = 1.3, at r = 0.18. Fig. |2] shows how the cor¬ 
responding diffusivities D depend on volume fraction 
and temperature. One sees that for T - 0.18, (p - 1.15 
gives essentially the lowest diffusivity; for </> = 1 one 
has standard behaviour, with D decreasing as density 
increases, while for (p - 13 the trend is reversed and 
one has anomalous behaviour. Fig. [T] suggests that the 
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dynamics at<p - I and at(p - 1.3 are very similar, with 
the two mean square displacement curves visually in¬ 
distinguishable; we return to this point later on. 

It is of particular interest to investigate whereas 
elastic models of the glass transition capture the ob¬ 
served density anomalies. Indeed, density anomalies are 
mainly driven by the density dependence of the entropy 
of the system, while elastic models are based on a purely 
energetic interpretation of the dynamics. 



<!> 


Figure 2: Volume fraction dependence of the isothermal dilfusivity, 
for dilferent values of the temperature as indicated. 


4. Inherent and liquid shear elasticity 

Elastic models III 0 of the glass transition are pred¬ 
icated on the assumption that the dynamics of super¬ 
cooled liquids consists of a series of jumps between 
local potential energy minima, i.e. inherent structures. 
The local curvature of these minima, which is related to 
the shear modulus of the system, sets the energy scale 
of these events and thus controls the relaxation dynam¬ 
ics. This energy-landscape interpretation makes it of 
interest to investigate the relation between the elastic 
properties measured in the liquid phase, and those ob¬ 
served after quenching the system to its inherent struc¬ 
ture. We have therefore measured the probability distri¬ 
bution function of the shear stress cr^y and of the in¬ 
herent shear stress, cr^y, whose variances are propor¬ 
tional to the instantaneous shear moduli. Fig. [3] shows 
that these distributions have a Gaussian-like shape at all 
temperatures. In systems interacting via inverse power- 
law potentials |@] (cr^y) decreases on cooling, and ap¬ 
proaches the value of ((cr^^)^) which is temperature in¬ 
dependent. We find that also in our system, (cr^y) de¬ 
creases on cooling. However, as is clear from Fig. 01 





Figure 3: (color online) Probability distribution of the sheai* stress and 
of the inherent shear stress, at ^ = 1.15, for dilferent temperatures as 
indicated. 



Figure 4: (color online) Temperature dependence of the variance of 
the shear stress for the parent liquid and the inherent stnactures, for 
^ = 1 (left panel), ^ = 1.15 (central panel) and (p = 1.3 (right panel). 


significant differences appear in the behaviour of the in¬ 
herent structure stresses, and the relative size of the in¬ 
stantaneous and inherent structure stresses. 

We note first that ((crj.p^) depends on the tempera¬ 
ture of the parent liquid. This behaviour signals the 
fact that the system explores different regions of the en¬ 
ergy landscape at different temperatures. In similar sys¬ 
tems, this feature has recently been exploited to prove 
that the volume fraction at which the inherent structure 
shear stress vanishes, the jamming volume fraction, de¬ 
pends on the temperature of the parent liquid d. This 
proves that the jamming volume fraction is protocol de¬ 
pendent IIITL IlST . At high temperatures the dynamics 
of the system is no longer affected by the energy land¬ 
scape, and ((cr^®)^) should be temperature independent. 
Fig. 0] suggests that we have reached this high tempera- 
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Figure 5: (color online) Temperature dependence of the average con¬ 
tact number per particle for the parent liquid (top) and the inherent 
structures (bottom), for (j> = 1.00. 


ture regime for ^ - 1. 

Looking next at the relation between instantaneous 
and inherent structure stresses, we find consistently (see 
Fig. nil that ((cr^)^) > {cP'^y), while the opposite relation 
is found in inverse-power law liquids. We can make 
sense of this result by considering that the fluctuations 
of the shear stress are related to the instantaneous shear 
modulus. For harmonic potentials, whose second (ra¬ 
dial) derivative of the potential is constant, the instanta¬ 
neous shear modulus is expected to scale with the den¬ 
sity of contacts, G oc. up where z is the average contact 
number per particle. Fig. |5] shows the temperature de¬ 
pendence of z for the parent liquid, where z increases on 
cooling, and for the inherent structures, where z a; 5.5 
is constant. At all temperatures, the average contact 
number of the inherent structures is larger than the av¬ 
erage contact number of the parent liquid. This implies 
that the inherent shear modulus 

larger than the shear shear modulus of the parent liquid. 
Accordingly, ((cr^®)^) > {crly), which 
is the relation we saw in Fig.|4] In power law potential 
systems the situation is different; the single particle bulk 
modulus, which is directly related to the second deriva¬ 
tive of the potential, is not constant but increases with¬ 
out bound as the interparticle distance becomes smaller. 
Thus, if the average interparticle distance increases due 
to more efficient packing when the system is quenched 
to its the inherent structure, then the shear modulus of 
the inherent structures is expected to be smaller than that 
of the parent liquid. 



Figure 6: (color online) Transient shear modulus and transient inher¬ 
ent structure shear modulus, for ^ = 1 and T = 0.16 (left panel) and 
T = 0.13 (right panel). 

5. Elastic models 

We now consider whereas the slowdown of the dy¬ 
namics is well described by elastic models of the glass 
transition. First we consider Dyre’s shoving model, ac¬ 
cording to which logT oc GpVatlksT, where as before 
Gp is the plateau shear modulus, and V^t a local ac¬ 
tivation volume. As usual we assume that the latter 
does not change with temperature, but now add the as¬ 
sumption that it is also independent of volume fraction. 
Fig. |6] shows the relaxation dynamics of the instanta¬ 
neous shear stress: the transient shear modulus G(f) ex¬ 
hibits the two-step decay typical of glasses. From G(f), 
we can estimate a plateau modulus Gp in the deeply su¬ 
percooled regime as the value of G(f) at the (intermedi¬ 
ate) time at which the derivative of G(f) with respect to 
log(f) is minimal. The stress relaxation time r can then 
be extracted via a stretched exponential fit of the final 
decay of G{t). 

Having obtained data for Gp and t, we can then as¬ 
sess the applicability of Dyre’s shoving model to our 
system: see Fig. |7] Here we show (black circles) that 
the relaxation time is correlated fairly tightly with the 
plateau shear modulus Gp In particular, the expected ex¬ 
ponential relation between relaxation time and plateau 
modulus divided by T is observed in the deeply super¬ 
cooled regime. It is important to note here that the figure 
combines points taken at different temperatures and vol¬ 
ume fractions. In particular, the filled black circles refer 
to r = 0.14 and volume fractions around the one giv¬ 
ing the minimal diffusivity. This proves that elasticity 
as measured by the plateau shear modulus is well corre¬ 
lated with the dynamics in the anomalous region. 

While the prediction of the shoving model appears to 
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Figure 7: (color online) Correlation of the relaxation time t with the 
plateau shear modulus Gp, the instantaneous shear modulus G(0), and 
instantaneous inherent structure shear modulus G*®. Data for a range 
of different temperatures and volume fractions are shown. The filled 
symbols refer to T = 0.14 and volume fractions (p in the region of the 
density anomaly. 


be reasonably well verified, its interpretation in terms of 
the features of the energy landscape of the system is not. 
Indeed, in this interpretation the plateau shear modulus 
should be related to the inherent structure shear mod¬ 
ulus, as recently found in power-law liquids 1113]. In 
these systems G'®(0) Gp and at all times G*®(f) < G(f) 
Intuitively, the instantaneous stress has larger fiuctua- 
tions over short time scales than the inherent structure 
stress, but around the timescale of the plateau the fast 
fiuctuations have averaged out and the relaxations of 
instantaneous and inherent structure stress track each 
other. 


In our short range harmonic repulsive systems, on 
the other hand, we find G*^(0) > G(0) as shown in 
Fig.|6] and as explained in the previous section. Since 
G(0) > Gp, this implies also G*®(0) > Gp'. the approxi¬ 
mate equality between these quantities no longer holds. 
In fact we find that G'®(0) is not even proportional to Gp. 
This is clear indirectly from Fig.|7] as no data collapse 
is found when t is plotted versus G^^{Q)/ksT instead 
of GpIk^T', the same holds if we use G(0) instead of 
G*®(0). This result clarifies that, while the dynamics of 
harmonic systems are determined by their elastic prop¬ 
erties, these properties are not related in a simple way to 
those of the energy landscape. This is consistent with a 
previous entropic interpretation of the observed density 
anomalies 1121] . 


We next consider, as a test of the second elastic model 
mentioned in the introduction, whereas the relaxation 
dynamics is also closely correlated with the DW factor. 


Figure 8: (color online) Time evolution of the exponent ^(0 chai'ac- 
terizing the diffusivity of the system, oc at different tem¬ 

peratures for volume fraction 0=1. 


Fig.m suggests that this may well be the case, as for 
T = 0.18 we saw that the mean square displacements 
at 0 = 1 and at 0 = 1.3 coincide to good accuracy at 
all times, consistent the macroscopic diffusivities also 
being equal between these two volume fractions. We 
have therefore investigated the existence of correlations 
between the diffusion coefficient D and the DW factor, 
{u^}. Operatively ||3], we define {u^} - {r~{tow))^ where 
tow is the time of minimal diffusivity. We determine 
this time by considering the time dependence of the dif¬ 
fusivity exponent b{t) - 51og((r“(f)))/51og(f), which 
varies in time from the value b -2, characteristic of the 
short time ballistic motion, to the value b - \ for the 
long time diffusive motion. In the supercooled regime, 
the mean square displacement at intermediate times is 
sub-diffusive, and b{t) < 1 has a minimum at some time 
tow- 

Fig. |9] displays the resulting dependence of the dif¬ 
fusion coefficient D on In the deeply supercooled 
regime of low D, the figure suggests that D is uniquely 
determined by (m^). At higher temperatures this is no 
longer the case. We note, however, that at higher tem¬ 
peratures the identification of (m^) is subject to large er¬ 
rors as the subdiffusive regime disappears. In addition, 
in the anomalous volume fraction range we observe the 
presence of a long subdiffusive regime with a nearly 
constant subdiffusive exponent, as illustrated in the inset 
of Fig. HI 

We conclude this section by returning to the point that 
the two elastic models we have considered are not in¬ 
dependent because {u^) is (mainly) determined by the 
shear elasticity Indeed, we do also find a clear 

correlation between the relaxation time and the DW fac- 
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6. Conclusions 



Figure 9: (color online) Parametric plot of the diffusion coefficient 
versus the Debye-Waller factor, for <^ = 1,1.15 and 1.3. In the deeply 
supercooled regime the diffusivity appears to be determined by (u^). 
At higher temperatures the determination of (u^) is affected by sub¬ 
stantial uncertainty because the subdiffusive regime disappears; in ad¬ 
dition, at high volume fraction {r^(t)) can develop an extended subd¬ 
iffusive regime as illustrated in the inset. 



Figure 10: a (color online) Con'elations between (a) relaxation time 
and DW factor, (b) DW factor and plateau modulus, and (c) diffusivity 
and relaxation time. The straight line in panel (c) is D oc with 
q K 0.72. Data points shown cover a range of different temperatures 
and volume fractions. 


tor, which is consistent with the Hall-Wolynes equation, 
T oc exp(a^/2(i/^)) as shown in Fig.fTOk. The combined 
validity of this relation and of Dyre’s model implies the 
relation T/Gp oc {u^), which is also compatible with 
our numerical data as Fig.fTOb shows. Finally, we note 
(Fig. [Toll that we also find a relation between diffusivity 
and relaxation time, with D oc and q x 0.72. 


We have demonstrated that elastic models of the glass 
transition correctly capture the slow dynamics of har¬ 
monic particle systems, in the temperature and volume 
fraction region where density anomalies are observed. 
However, we have found the relevant elastic constants 
not to be related to the features of the energy landscape 
of the system. This result challenges the usual poten¬ 
tial energy landscape interpretation of elastic models, 
and suggests that they should generally be interpreted 
by referring to the free energy landscape of the system. 
Future directions include an investigation of the valid¬ 
ity of elastic models in other liquids with density driven 
anomalies, such as the Gaussian potential, the Jagla po¬ 
tential, or water-like model, as well as in liquids with 
temperature driven anomalies such as the sticky hard- 
sphere model. 
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